!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
recursive subroutine hilbert_curve (parent, current, parent_cell, pos, order, location, var_regrid_flag)
!==================================================================================================================================|
!  order: Number of iterations
!  location: The location of current block relative to the left, bottom, behind corner in current hilbert order

   use parameters
   use block_define
   implicit none

   integer(4), intent(in) :: parent_cell(n_dims, n_cell), pos, order, var_regrid_flag, location(n_dims)
   type(block), pointer, intent(in) :: parent
   type(block), pointer, intent(inout) :: current

   integer(4) :: i_cell, i_dim
   integer(4) :: entry_exit_cell(n_dims, 2, n_cell), gene_cell(n_dims, 2), unit_cell_tmp(n_dims, n_cell)
   integer(4) :: location_tmp(n_dims), mirror(n_dims), exchange_pos1, exchange_pos2

!----------------------------------------------------------------------------------------------------------------------------------|
!  confirm the entry and exit for all 8 cells
!----------------------------------------------------------------------------------------------------------------------------------|
   do i_cell = 1, n_cell
      if (i_cell .le. n_cell / 2) then
         if (i_cell .eq. 1) then
            entry_exit_cell(:, 1, i_cell) = parent_cell(:, i_cell)
            entry_exit_cell(:, 2, i_cell) = entry_exit_cell(:, 1, i_cell) .neqv.                                                   &
               (parent_cell(:, i_cell) .neqv. parent_cell(:, i_cell + 1))
         else
            entry_exit_cell(:, 1, i_cell) = entry_exit_cell(:, 2, i_cell - 1) .neqv.                                               &
               (parent_cell(:, i_cell - 1) .neqv. parent_cell(:, i_cell))
            if (count ((entry_exit_cell(:, 1, i_cell) .and. (parent_cell(:, i_cell) .neqv. parent_cell(:, i_cell + 1))) .eq.       &
                     (parent_cell(:, i_cell) .and. (parent_cell(:, i_cell) .neqv. parent_cell(:, i_cell + 1)))) .eq. n_dims) then
               entry_exit_cell(:, 2, i_cell) = entry_exit_cell(:, 1, i_cell) .neqv.                                                &
                  (parent_cell(:, i_cell) .neqv. parent_cell(:, i_cell + 1))
            else
               entry_exit_cell(:, 2, i_cell) = entry_exit_cell(:, 1, i_cell) .neqv. cshift ((parent_cell(:, i_cell) .neqv.         &
                  parent_cell(:, i_cell + 1)), 1)
            endif
         endif
      else
         mirror = parent_cell(:, n_cell / 2) .neqv. parent_cell(:, n_cell / 2 + 1)
         entry_exit_cell(:, 1, i_cell) = entry_exit_cell(:, 2, n_cell - i_cell + 1) .neqv. mirror
         entry_exit_cell(:, 2, i_cell) = entry_exit_cell(:, 1, n_cell - i_cell + 1) .neqv. mirror
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  generate the gene sequence for the target cell (defined by pos)
!----------------------------------------------------------------------------------------------------------------------------------|
   gene_cell(:, 1) = (unit_cell(:, 1) .neqv. unit_cell(:, n_cell)) .neqv.                                                          &
                        (entry_exit_cell(:, 1, pos) .neqv. entry_exit_cell(:, 2, pos))
   gene_cell(:, 2) = unit_cell(:, 1) .neqv. entry_exit_cell(:, 1, pos)

!----------------------------------------------------------------------------------------------------------------------------------|
!  exchange the unit cell to the new cell according to the gene sequence (this operation is executed only once)
!----------------------------------------------------------------------------------------------------------------------------------|
   unit_cell_tmp = unit_cell
   exchange_pos1 = -1
   do i_dim = 1, n_dims
      if (gene_cell(i_dim, 1) .eq. 1) exchange_pos1 = i_dim
   enddo

   exchange_pos2 = -1
   do i_dim = n_dims, 1, -1
      if (gene_cell(i_dim, 1) .eq. 1) exchange_pos2 = i_dim
   enddo

   if (exchange_pos1 * exchange_pos2 .lt. 0) then
      write(*, *) ' Can not execute the exchanging operation.'
      call mpi_abort (mpi_comm_world, 1, error)
   endif

   if (exchange_pos1 .ne. -1 .and. exchange_pos2 .ne. -1) then
      unit_cell_tmp(exchange_pos1, :) = unit_cell(exchange_pos2, :)
      unit_cell_tmp(exchange_pos2, :) = unit_cell(exchange_pos1, :)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  reverse the unit cell to the new cell according to the gene sequence (this operation can be executed more than once)
!----------------------------------------------------------------------------------------------------------------------------------|
   do i_dim = 1, n_dims
      if (gene_cell(i_dim, 2) .eq. 1) then
         unit_cell_tmp(i_dim, :) = 1 - unit_cell_tmp(i_dim, :)
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  Iterate the next generation or creat the links if iteration stopped
!----------------------------------------------------------------------------------------------------------------------------------|
   if (order .gt. 1) then
      do i_cell = 1, n_cell
         call hilbert_curve (parent, current, unit_cell_tmp, i_cell, order - 1,                                                    &
            location + unit_cell_tmp(:, i_cell) * 2 ** (order - 1), var_regrid_flag)
      enddo
   else
      do i_cell = 1, n_cell
         location_tmp = unit_cell_tmp(:, i_cell) + location
!----------------------------------------------------------------------------------------------------------------------------------|
!        creat the links for the block in the computational box
!----------------------------------------------------------------------------------------------------------------------------------|
         if (location_tmp(1) .le. ini_block_x .and. location_tmp(2) .le. ini_block_y .and. location_tmp(3) .le. ini_block_z) then
            if (var_regrid_flag .eq. cpu_id) then
               current%next => parent%child(location_tmp(1), location_tmp(2), location_tmp(3))%pblock
               current%next%last => current
            endif
            current%framework_next => parent%child(location_tmp(1), location_tmp(2), location_tmp(3))%pblock
            current%framework_next%framework_last => current
            current => current%framework_next
            current%parent_cell = parent_cell
            current%hilbert_pos = i_cell
         endif
      enddo
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine hilbert_curve
